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We present a statistical field theory to describe large length scale effects induced by solutes in a 
cold and otherwise placid liquid. The theory divides space into a cubic grid of cells. The side length 
of each cell is of the order of the bulk correlation length of the bulk liquid. Large length scale 
states of the cells are specified with an Ising variable. Finer length scale effects are described with 
a Gaussian field, with mean and variance affected by both the large length scale field and by the 
constraints imposed by solutes. In the absence of solutes and corresponding constraints, integration 
over the Gaussian field yields an effective lattice gas Hamiltonian for the large length scale field. In 
the presence of solutes, the integration adds additional terms to this Hamiltonian. We identify these 
terms analytically. They can provoke large length scale effects, such as the formation of interfaces 
and depletion layers. We apply our theory to compute the reversible work to form a bubble in 
liquid water, as a function of the bubble radius. Comparison with molecular simulation results 
for the same function indicates that the theory is reasonably accurate. Importantly, simulating 
the large length scale field involves binary arithmetic only. It thus provides a computationally 
convenient scheme to incorporate explicit solvent dynamics and structure in simulation studies of 
large molecular assemblies. 

PACS numbers: 61.20.-p, 61.20.Gy, 68.08.-p, 82.70.Uv, 87.15.Aa 



I. INTRODUCTION 

We have constructed a tractable model for describing 
density fluctuations in a cold liquid at both small and 
large length scales. The model allows us to analyze at 
a microscopic level the effects of solvated surfaces and 
large molecular assemblies, perhaps of biophysical rele- 
vance. This paper presents the model and demonstrates 
its tractability. 

A cold liquid is a fluid that is well below the criti- 
cal temperature. Water at ambient conditions is an ex- 
ample. When unperturbed, it will have no significant 
large length scale fluctuations. It is nearly incompress- 
ible. When perturbed by a sufficiently extended surface, 
however, a cold liquid may exhibit large length scale fluc- 
tuations, akin to a phase transition, in the vicinity of the 
surface. This phenomenon occurs when another phase 
is close to coexistence with the liquid, and when inter- 
actions with the surface favors the other phase over the 
liquid. This coincidence of conditions is pertinent, for 
instance, to hydrophobic effects. In particular, water 
at ambient conditions lies close to coexistence with its 
vapor. Further, the dcmixing of oil and water and the 
associated large oil-water surface tension indicates that 
a large hydrophobic (i.e., oily) surface favors vapor over 
liquid water. 

Indeed, Lum, Chandler and Weeks (LCW ) ||] have 
demonstrated that oily surfaces extending over 1 nm 
or more will nucleate a layer of depleted water density 
and concomitant large length scale correlations. In con- 
trast, perturbations from smaller hydrophobic surfaces, 
less than 1 nm across, do not nucleate such a drying layer 
and affect only small length scale fluctuations in the liq- 
uid. Since hydrophobicity vividly manifests the interplay 
and competition between small and large length scale 



fluctuations in a cold liquid, we have chosen in this pa- 
per to focus attention on it. One bencflt of our analysis 
is an understanding of the results of LCW theory from 
a perspective that is numerically simpler and physically 
more transparent than the original LCW development. 
Generalizations of our approach to other phenomena, in- 
cluding the effects of strong associative interactions be- 
tween solutes and solvent, should be apparent. 

The main idea of our approach is to create a statistical 
field theory where the molecular density fleld is decom- 
posed into two parts. One part varies on large length 
scales only. The other varies on small length scales. For 
a cold fluid that is homogeneous and therefore nearly 
incompressible, the large length scale field is nearly con- 
stant and equal to the mean density of the bulk liquid. 
Even for this homogeneous case, however, small length 
scale fluctuations are always present. To a remarkable 
extent |H,D, the statistics of these fluctuations is Gaus- 
sian with a variance determined by the structure factor of 
the bulk liquid. Accurate molecular theories of solvation 
and liquid structure at small length scales - the Percus- 
Yevick equation for hard sphere fluids the mean 

spherical approximation , the Pratt-Chandler theory 
of hydrophobicity |^ , the reference interaction site model 
(RISM) |^,|o| - are consequences of such statistics 0. 
This Gaussian statistics for small length scale fluctua- 
tions is an important element of the weight functional 
(or Hamiltonian) we construct. These fluctuations are 
coupled, of course, to the large length scale density field, 
and they are also constrained by the presence of solutes. 
Due to the coupling and constraints, the variance of the 
small length scale fluctuations may differ markedly from 
that of the homogeneous bulk fluid. 

The Hamiltonian for our model is presented in Sec. ||. 
The large length scale density fleld supports possible 
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phase coexistence and interfaces. As such, we see in that 
section how the coupUng between small and large length 
scale fields may lead to solute induced interfaces in a cold 
fluid. Our treatment of this coupling is inspired by the 
work of Weeks and his coworkers y,f2[. They related 
the coupling to unbalanced attractive forces that result 
from local inhomogeneities in the fluid. Analytical inte- 
gration over the small length scale field is possible due to 
its Gaussian statistics. The integration yields an effective 
Hamiltonian functional for the large length scale field. In 
Sec. Ill, we describe how this integration can be used to 
study solvation. This step also lays the foundation for a 
numerical scheme where the solvent is simulated at the 
level of the large length scale density field. Such a scheme 
involves only binary arithmetic and is much more efficient 
than an atomic level simulation. In fact, it is sufficiently 
efficient to make possible the study of phenomena like 
self-assembly of biological structures. 

In Sec. 1^, we discuss the results of our treatment in 
various limits. In the absence of any solutes, the effec- 
tive Hamiltonian for the large length scale density corre- 
sponds to the lattice gas model |jl^. In the presence of 
solutes that are small in size and number, only those den- 
sity fiuctuations at small length scales are relevant, and 
our model reduces to the Gaussian model of Pratt and 
Chandler and the closely related information theory 
approach of Hummer, Pratt and cowor kers 1]. In the 
presence of large solutes, a mean field approximation to 
our model coincides with the LCW theory fy]. 

A numerical application is given in Sec. We first 
show how the parameters in our model can be estimated 
from experimentally accessible quantities. We then ex- 
plicitly treat the solvation of an ideal hydrophobic sphere 
in water and compare our results with those of an atom- 
istic simulation iQ. Finally, in Sec. VI, we discuss im- 



plications and possible extensions of this work. 



II. MODEL 

Figure |l| illustrates the essential features of a cold fluid 
in the presence of a solute. The solute is of arbitrary 
size and shape. If it is small, the solvent will wet its 
surface. In contrast, if the solute is large with extended 
hydrophobic surfaces, solvent density near the solute will 
be depleted relative to the density of the bulk liquid, pi. 
This drying-like phenomenon occurs because the solvent 
experiences significant unbalanced attractive forces near 
the hydrophobic surface. These forces induce depletion. 
The solute can also have patches of associative interac- 
tions. Adjacent to those patches, the molecular density 
of the solvent will be close to or perhaps greater than 
that of the bulk liquid. 

In our description of solvation, we make a distinc- 
tion between strong forces and weak solvent-solute forces. 
The repulsive nearly hard core interactions between so- 
lute and solvent molecules are strong forces. So too are 
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p(r)=pi+5p(r) 



FIG. 1. Sketch of our model of a cold liquid in the pres- 
ence of a solute. The solute excludes a volume v^^ from the 
solvent (black region) and has a hydrophilic patch of strong 
associative interactions with the solvent; the patch imposes 
a constraint on the solvent density in the volume Vi, (gray 
area). The solvent is divided into cells of width I; each cell 
is either filled with liquid {rii — 1) or vapor {ui = 0). The 
field rii describes density fluctuations on length scales larger 
than the lattice spacing. This field supports phase transitions. 
Density fluctuations on length scales smaller than the lattice 
spacing are described by the field Sp{r). This field describes 
molecular detail such as the highly oscillatory profiles for the 
average density near small solutes. We thus write the density 
as p{ri) — Uipi + Sp{ri). The field 5p{r) is assumed to obey 
Gaussian statistics. 



associative interactions between solute and solvent. On 
the other hand, dispersion interactions between solute 
and solvent molecules are weak forces. In some cases, 
electrostatic forces are weak forces. Weak interactions 
are described in our treatment by an interaction poten- 
tial acting between the solute and the solvent density. 
In contrast, strong forces are treated according to the 
constraints they impose upon the solvent density fluctu- 
ations. For example, the effect of a solute repulsive core 
is mainly to exclude solvent from a volume, Vex m Figj^. 
The effect of these forces may be described as a constraint 
permitting only those fluctuations in the solvent density, 
p(r), that leave Vcx empty of solvent, i.e., p(r) = for 
r G I'cx 0] ■ Similarly, associative interactions may cause 
n water molecules to be bound within a specific region 
close to the solute. In Fig.|l|, this region is Ub. This effect 
can be treated by constraining the integral of p(r) over 
the volume vt, to equal n We do not apply this lat- 
ter idea in the current paper, although the methods by 
which it can be implemented should be clear from our 
treatment of the former. 

Since the fluid is assumed to be cold, the regions of gas 
or vapor can be clearly distinguished from those of liquid. 
The density of the vapor is typically orders of magnitude 
smaller than that of the liquid. In such a situation, it is 
natural to divide space into a grid of cells, where each 
cell contains either gas or liquid. We use cubic cells, and 
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take the distance across each cell, I, to be comparable 
to the bulk liquid correlation length, ^. In that case, a 
binary choice of states within a cell, either gas or liquid, 
provides a reasonable coarse grained rendering of likely 
configurations of the fluid. We can thus define a field, rii, 
that takes on the value of 1 if cell i contains liquid and 
if it contains gas. The molecular density we associate 
with this field is riipi. This field, n^, or equivalently riipi, 
is the large length scale field in our model. It can be 
used together with a second field, Sp{r), to describe the 
density on length scales both larger and smaller than I. 
In particular, for positions r within cell i, r^, we write 
the net density as 

p{ri) = Uipi + 5p{ri). (1) 

All of space is spanned by the set of r^, i.e., / dr = 

While the field nj is binary and can be used to de- 
scribe a liquid-gas phase transition, the field Sp{r) has a 
very different character. It supports neither phase tran- 
sitions nor interfaces, but it does describe small length 
scale structures such as those manifesting the granularity 
of the solvent (e.g., the oscillatory profiles of the average 
liquid density in the vicinity of a small solute) . It must 
be possible, therefore, that 6p{r) can take on a variety of 
values. As indicated in the Introduction, it is a reason- 
able approximation to adopt the simplest possible statis- 
tics for this field. Namely, we assume it is Gaussian and 
define its variance to be 

x[r.,r;.;{n4] = (Mr.)'5p(r;-)>{n.}. (2) 

Here, (■■•){„j.} indicates the ensemble average over den- 
sity fiuctuations for a given configuration of the field nk- 
The dependence upon Uk is significant. When Uk — 1 for 
all k, corresponding to a cold liquid with absolutely no 
large length scale fluctuations, Sp{ri) has zero mean, and 
its variance reduces to the response function of the bulk 
fluid, 

X{ri,r'j;pi) = pi5{r, - r'^) + pfh{\r, -r^|;p/), (3) 

where /i(|ri — r^- 1; p/) -I- 1 is the radial distribution func- 
tion of the uniform fluid at density pi. On the other hand, 
within a cell that contains vapor (rt,; = 0), small length 
scale fluctuations are very small. Our model employs the 
approximation that Spiji) = whenever rii = 0. There- 
fore, we imagine that Sp(r) is a Gaussian field, with a 
weight functional being that of the bulk fluid, but con- 
strained to be zero whenever Ui = 0. The response func- 
tion for such a field is 

X[ri,r^-; {uk}] = x(r.j, r^-; p,) 

k I •' 

xx^'[r'irT;{nk}]x{vTy,:Pi)- (4) 



Here, 

Xg[ri,r^; {nk}\ = x(r,,r^;p/) if n.i = n.j = 0, 

= otherwise, (5) 

is the (vijV'j) element of the matrix Xg- Similarly, the 
matrix with elements x^T^ [r*, r^; {n^}] is also non-zero 
only when Ui = Uj = 0. In that space, where 71^ = 0, 
Xg^ is the inverse of Xg- These relations project the 
matrix with elements x~^[^ii^'j! {^k}] onto the space of 
cells for which = 1. We adopt these relations to define 
our model of Gaussian statistics for Sp{r). 

With the lattice spacing as large as the bulk correla- 
tion length, the field rii pi is nearly incompressible. This 
means that, in the absence of any strong perturbations 
on the fluid, the fleld rii is essentially constant. In this 
case of the unperturbed (i.e., uniform) fluid, density 
fluctuations are described almost entirely by the field 
Sp{r). The compressibility of the uniform fluid is con- 
tained in the variance for Sp{r). In this context, con- 
sider the behavior of the bulk liquid structure factor, 
/ c?(r— r')x(r, r'; p/) exp[ik-(r— r')]. The long- wave length 
limit of the structure factor is proportional to the bulk 
compressibility. It approaches this limit with a plateau. 
In particular, for k values smaller than some finite wave- 
vector kc , the structure factor is essentially constant. 
The grid spacing we use to define large length scales co- 
incides with I 27r/fcc. 

To within a physically irrelevant metric factor, the par- 
tition function for our model is 

^ = E / T^Sp{r)C[{n,,},Sp{v)] 

X exp{~m{nk},5p{r)]), (6) 

where / V5p{r) = J IliT>Sp{ri) denotes the functional in- 
tegration over the small length scale field, H[{nk} , Sp{r)] 
is the Hamiltonian as a functional of both rii and Sp{r), 
and is Boltzmann's constant time temperature, k^T. 
The quantity C[{nk} , <^p(r)] is a constraint functional. 
It has unit weight when the field Sp{r) together with 
{rii} satisfy whatever constraints are imposed by strong 
forces, and it is zero otherwise. Since {rii} and Sp{r) have 
greatly different character, the summation and integra- 
tion in Eq. do not redundantly count configuration 
space to any significant degree. 

In our model, there are three principal contributions 
to the Hamiltonian, H[{nk} , Sp{r)]. One is a lattice gas 
Hamiltonian for the large length scale field, 

H-L[{nk}] = -fJ-^rii - e ^ riiUj. (7) 

i <i,j> 

Here, p is the imposed chemical potential, the sum la- 
beled with (ij) is over all nearest neighbor pairs of cells, 
and the interaction parameter e determines the energetic 
cost of creating a vapor-liquid interface. Importantly, the 
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lattice gas model supports phase transitions and sustains 
gas-liquid interfaces. 

A second contribution to the Hamiltonian ensures the 
Gaussian weight for the small length scale field. From 
the principle of equipartition, this contribution must be 



^ y dr, y dr;.<5p(r,)x-Mr.,r; ;K}]Mr;). (8) 



A third contribution to the Hamiltonian gives the cou- 
pling between the rii and Sp(r) fields arising from unbal- 
anced forces. According to the arguments provided by 
Weeks and coworkers for simple fluids ||l|,|l^, the unbal- 
ancing potential acting on for a simple fluid is well 
estimated by —2a{dp{r)) . Here, apf is the energy den- 
sity of the uniform liquid at density p;, and the overbar 
denotes a coarse-graining of the density fluctuation Sp{r) 
over a length scale comparable to the bulk correlation 
length. Based upon this estimate, we write the contribu- 
tion to the Hamiltonian from unbalanced forces as 



-2e' 



E 

<i,3> 



dvi 5p{vi) 



(9) 



For simple fluids with only one energy and length scale, 
e' = e. This equality implies that the lattice gas param- 
eters e and I are sufflcient to determine both the surface 
tension and the energy density of the liquid. For more 
complex fluids, including water, the effects of orienta- 
tional degrees of freedom may introduce multiple micro- 
scopic length scales, and as a result, e' could differ from 
e. This possibility was ignored in Ref. M, but will be 



examined in Sec. VI 



By combining all three contributions, we arrive at our 
final result for the Hamiltonian of our model. It is 



H[{nk},Spir)] 



H^[{nk}]~2e'J2 J dnSp{r. 
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<l:J> 



Pll"" 



i?„o™[{nfc}], (10) 



where 



X (j)k x[ri,r';{nfc}] (f>i 



and 



fceTln Vdd^, 



PiP 



(11) 



(12) 



Here, the quantity detx is the determinant of the ma- 
trix with elements r^; {"fife}]. The last term in 
Eq. (p^, -ff„orm[{?T-fc}], provides a normalization constant 
for the functional integration over 5p{r). When there 
are no strong forces, so that the constraint functional 
C[{nfc} , Sp{r)] is simply unity, the effective Hamiltonian 
for the Hi field should be exactly the lattice gas Hamil- 
tonian, Eq. (0). The last term in Eq. ( p^ ensures that 
the integration over Sp{r) for this case will indeed yield 
this result. 



III. THEORY OF SOLVATION 

The excess chemical potential of a solute, Ap, is given 
by @ 

(3Ap = -In^ = -ln(exp(-/?f/s))o. (13) 

Here S is the partition function for the unperturbed sol- 
vent and Ss is the partition function for the system in 
the presence of a (fixed) solute. The energy Us is the 
energy of interaction between the solute and the solvent 
molecules and the subscript denotes an ensemble aver- 
age over the unperturbed solvent. 

For simplicity, we will consider the solvation of an ideal 
hydrophobic solute in water - a particle that excludes wa- 
ter from a region Wox, but has no other interactions with 
the solvent. A hard sphere is an example of an ideal 
hydrophobic solute. It excludes solvent from a volume 
= (4/3)7ri?'^, where the radius R is the distance of 
closest approach between water and solute. The parti- 
tion function of the system in the presence of such a so- 
lute is equal to the partition function of the unperturbed 
solvent, but with the constraint that no solvent exists in- 
side the excluded volume. In other words, the constraint 
functional for this case is 



C[{nk},6p{r)]^ Yl S[nipi+6p{ri)]. 



(14) 



Accordingly, the partition function in the presence of an 
ideal hydrophobic solute is 



E / U^^pi^^ 




d[nipi + dp{r.i)] 



X exp{~/3H[{nk},Sp{r)]). 



(15) 



For an ideal hydrophobic solute, the ratio of partition 
functions Ss/S, equals the probability of observing no 
solvent molecules inside the volume v^^. Equivalently, 
it corresponds to the probability of observing a cavity 
of volume v^^ inside the solvent; it is also equal to the 
probability that a solute can be inserted into the solvent 
without creating any overlap with the solvent molecules. 
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The excess chemical potential of an ideal solute could be 
obtained by imposing an alternative constraint, 



drpivi 



(16) 



Were our treatment completely consistent with the par- 
ticulate nature of matter, the two constraints, as given 
by Eq. (|lj) and Eq. (p^), would be equivalent. But in 
fact, Gaussian statistics for the field (5p(r) cannot be com- 
pletely consistent with this nature of matter, and the 
two constraint functionals will yield somewhat different 
results. 



To evaluate the partition function, Eq. (15), it is conve- 
nient to rewrite the constraint functional with the Fourier 
representation of delta functions. Namely, 



^s = E /n^'5p(r,) jY[V^{v,) 



exp -PH[{nk},6p{r)] 



driV(ri) [uipi + Sp{ri)] 



(17) 



Functional integration over both Sp{r) and ^p{r) is now 
straightforward, yielding 



(18) 



where the effective Hamiltonian H[{nk}] is 
H[{nk}] ^ Hi^[{nk}] 

X xrJ-i^u^'j-Ank}] {njpi + /(r'j)) 

(19) 



fceTln Vdetx.„, 



with 



fc(nnj) 



"-fc - 1 
Pll^ 



X[r,,r';{nfe}]. (20) 



Here, Xi^^ has elements x[ri, r^; {rifc}] for and both 
within the excluded volume, and a sum over fc(nnj') is 
over cells k that are nearest neighbors to cell j. 

The evaluation of i?[{rife}] requires the calculation of 
various integrals and matrix inverses. These quantities 
can be conveniently estimated to a good approximation 
by exploiting the fact that the lattice spacing is on the or- 
der of the bulk correlation length. In particular, since the 
bulk correlation function, x(rz, r^-; Pi) vanishes quickly for 
jr — v'\ larger than that length, x[ri, ; {rik}] as given by 
Eq. (Q), can be approximated by 



x[Y^,r'^■,{nk}] = xiri,r'j;pi), for 
= 0, otherwise. 



(21) 



Furthermore, the relatively large size of the cells allows 
us to restrict the sum in Eq. (GG) to the i = j term, and 
to take the integral over all space, rather than over one 
cell. As such, we arrive at a much simplified form for 
/(rj, and therefore. 



/ Pi 

niVi€ K — 



E(^fe-l)- (22) 

k{nni) 



Here k is the isothermal compressibility of the uniform 
fluid, which is related to the response function via k = 
Px{0)/p'i, where x(0) is the long- wave length limit of the 
Fourier transform of the structure factor. Note that fi is 
zero, when cell i is not liquid, i.e., when rii = 0. 

Finally, with x(ri,r^;/0/) provided as input into the 
theory, we must choose a set of basis functions that span 
the space of the excluded manifold. This allows us to 
perform the inversion of Xin[ri, r^-; {rifc}] in the represen- 
tation prescribed by that basis. We use the approxima- 
tion of one basis function spanning the excluded volume 
and to take Xi„[ri, r^; = xi^ii'^'j', Pi) for all cells i 

inside the excluded volume We then arrive at our 

principal result: 



H[{nk}] = H^[{nk}] 
+ kBT 

ij(occ) 



rii [piVi + fi] [piVj + fj] Uj 



2ct„<, 



A:bT In •\/27r(T^; 



Hi^[{nk}]+Hs[v^^;{nk}]. 



(23) 



Here the sum over i, j(occ) is over over cells i and j that 
are occupied by the solute; Vi is the volume occupied by 
the solute in cell i, and 



dr dr'x(r,r';p/). 



(24) 



In the one-basis set approximation for xl^i, r^-; {"■*:}], em- 
ployed to arrive at Eq. (^3|) , the effect of the constraint 
functional as given by Eq. (^^, reduces to that of the 
constraint functional as given by Eq. (p^. 

The term Hs[v^^;{nk}] contains all the effects of the 
interaction between the solute and the ideal hydrophobic 
solvent. It increases with increasing solute size, if — 1 
for the cells i that are occupied by the solute. The in- 
teraction term solely arises from the constraint that is 
imposed on the allowed density fluctuations of the sol- 
vent. This idea, that solvation of a hydrophobic species 
is equivalent to the effect of imposing a constraint on 
the solvent density, is an important feature of our model. 
Interestingly, the excess chemical potential of the solute 
can be obtained by averaging this interaction free energy 
as follows: 
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In- 



E{„,|exp(-/3g[{nfc}]) 
S{«.}exp(-/3i?L[{nfc}]) 
= -ln(cxp(-/3i/s[fox;{nfc}])L, 



(25) 
the 



where indicates the ensemble average with 

Hamiltonian Hi,[{nk}] = H[{nk}] - -ffs Kgii"-fe }] ■ 

The simple formula for Eq. (P3[), and simi- 

lar formulas for more general cases, can be of enormous 
practical benefit for studying self-assembly. Such stud- 
ies usually require large system sizes. In those cases, the 
treatment of solvent is a primary computational bottle- 
neck. This is because large solutes are solvated by huge 
number of solvent molecules, and an atomistic treatment 
involves a correspondingly large number of coordinates 
and momenta. The formula for H[{nk}], however, lays 
the foundation for a scheme in which only the solutes are 
treated explicitly at the atomic level; the solutes can be 
moved by a continuous Monte Carlo or molecular dynam- 
ics scheme. The solvent, on the other hand, is simulated 
in terms of the large length scale density field, rii. That 
field can be propagated by a dynamic Monte Carlo pro- 
cedure, manipulating only binary numbers. More details 
of this scheme will be discussed in a forthcoming publi- 
cation pTf . 



IV. LIMITING RESULTS AND COMPARISONS 
WITH OTHER THEORIES 

Consider first the case where rii — 1 for all i. This 
case is physically pertinent for solutes small in size and 
in number because the concomitantly small value of Vcx 
leads to relatively small free energetic costs for having 
rii — \ for all cells i, even for the cells that are occupied 
by the solute. Specifically, when the solutes occupy rela- 
tively small volumes, the amount that iJgiWox; {'^fc}] will 
decrease by changing Ui from 1 to will not compensate 
the corresponding increase in iJLll'T-fc}]- With = 1 for 
all i, Hi^[n] and i?„orm[{fii}] become constants and thus 
irrelevant. The response function xfr^, ; {rifc}] reduces 
to the response function of the uniform fluid r'; Pi)- 
Further, the coupling term in Eq. becomes identi- 
cally zero. As such, the Hamiltonian for the model re- 
duces to that of the Gaussian model of Pratt and Chan- 
dler d 11|, namely H [{rii} , <5p(r)] Hq [6p{v)] , where 



knT 



dr 



dr'Sp{r)x-\r,r';pi)Sp{r'), 
(26) 



with 5p{r) = p(r) — pi, and the response function 
X~^(r, r';p;) being the response function of the uniform 
fluid. Similarly, applying Eq. (p5|), we obtain the ex- 
cess chemical potential for the ideal hydrophobic solute: 



f3Ap{v^^) = -ln(exp(-/3i?sK.;{l}])L = PHs[v^^; {1}] 
But, if rii = 1 for all i, then fi — for all i, and so 



/3A/i(z;„) ~ p^i/2a„„^ + In y^2T^^, (27) 

For excluded volumes not unphysically small, this for- 
mula for the excess chemical potential is the solvation 
energy result of Hummer and Pratt and their coworkers 
We see that it is the result of Gaussian statistics 
in the one basis set approximation. 

In contrast to setting rii = \ for all i, the treatment of 
LCW 1^ assumes that rii can be replaced by {rii) , where 
(rii) is an estimate of the mean value of rii. In that case, 
Eq. (p5|) gives 



PAp{v,^ 



Hi^[{{nk)}] - Hi^m] + /3i?sbc.; {{rik)}] 
H^[{{n,)}]-H^[{l}] 



UbT ^ (n^) (rij) p^ViVj/2ay^ 

i,j{occ) 



(28) 



where the second approximate equality follows from 
the flrst after neglecting fi and fj in comparison with 
Vipi and Vjpi, respectively. Except in the crossover 
regime, this is usually a reasonable approximation, be- 
cause Kpi/P ^ 10~^. Within notational differences, 
Eq. (^8|) is the solvation energy result given by LCW 
theory, Eq. (9) of Ref. ^. 

The LCW formula for the mean large length scale field, 
(rii), can also be understood from our model. In particu- 
lar, in Eq. ( [lo| ) , let us replace in the second term the field 
Sp{ri) with its mean, (6p{ri)). With this replacement, 
the first two terms in Eq. (|lO| ) give the mean molecular 
field acting on nj: 



i(nnj) 

-p - e 



e(n,)+e' / dr, {Sp{r)) / pif 



J2 [{n..) + {Sp{v,))/pi 



(29) 



i(nnj) 



where the approximate equality follows principally from 
approximating e' with e. For the coarse graining indi- 
cated by the over-bar, we use I as the coarse graining 
length. Other contributions to the mean molecular field 
come from the quadratic term in Sp{r) and from H^„^. 
These, however, are small outside the crossover regime, 
either because they appear in the logarithm or because 
they arise from unlikely configurations, where one neigh- 
bouring cell is filled while another is empty. With the 
molecular field in Eq. ( p9| ) , the LCW self-consistent equa- 
tion for (uj) is obtained. Specifically, since both {rii) and 
{5p{vi))/pi vary slowly in space, they may be expanded 
for i close to j about {rij) and {5p{rj))/pi, respectively. 
Truncating the expansion of (rii) at the square gradient 
order, and the expansion for {5p{vi)) / pi at lowest order. 
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Eq. (H) gives Eq. (5) of Ref. §. Thus, the principal 
results of LCW can be understood as a mean field ap- 
proximations to the model we have presented herein. 



V. APPLICATION TO HYDROPHOBIC 
SOLVATION: FREE ENERGY TO HYDRATE A 
SPHERICAL BUBBLE 

A. Parameters 

To apply our model, experimentally accessible quanti- 
ties such as the structure factor and the surface tension 
of the solvent must be known. For the structure fac- 
tor, we use the data of Narten and Levy for water [ p^ . 
For the energy parameter e, we consider its connection 
to the experimental vapor-liquid surface tension of wa- 
ter, 7. Namely, since the liquid is cold, we use the low 
temperature relation. 



7 



2P' 



(30) 



With 7 = 70 mN/m, this yields e = 6.02 ksT. 

For the energy parameter e' , we consider its approx- 
imate connection with the energy density of the fluid, 
—apf. In particular. 
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FIG. 2. The excess chemical potential per unit area of a 
hard sphere in water as a function of its size; the hard sphere 
excludes water from a spherical volume of radius R. The dot- 
ted line indicates the average potential energy from the bare 
lattice gas model as given by //l [n] in Eq. (ES) and the dashed 
line gives the average potential energy of the solute-solvent in- 
teraction as given by Hs[n] in Eq. (^^. It is seen that the 
lattice artifacts in the two energy contributions tend to cancel 
each other. The horizontal line lies at the value of the surface 
tension 7 of the vapor-liquid interface of water. 



ze 



(31) 



where z = 6 is the coordination number of the lattice of 
the cubic lattice gas model |Q . Since the lattice spacing 
I should be on the order of the bulk correlation length, 
^, we have chosen ^ k, I — A.2 K. We have verified that 
the results, such as those given below, do not depend 
strongly on the precise value of the lattice spacing by 
varying it between Z = 3.5 A and I = 5.0 A. The energy 
density parameter a was derived from the internal energy 
U , which was obtained in the following way: 



U{2m K) « Ai?,,p(373 K) 



298 



CpdT, 



(32) 



373 



where AH^^^, is the heat of vaporization and Cp is the 
heat capacity. With AH^.^^ — —40.7 kJ/mol, Cp — 
75.3 J/K/mol, this yields a = 566 A^ and thus 
e' = 15.2 ksT. Note that e' > e. The inequality indi- 
cates that the unbalancing potential is strong enough to 
induce a vapor layer near the surface of a large hydropho- 
bic object. In fact, the results of the numerical applica- 
tion given below, suggest that the unbalancing potential 
is even stronger. 

For the imposed chemical potential, /i, we use 



^J■ = Mcoox + Afi w /i^ 



(33) 



where /^coox — ~3e is the chemical potential at coexis- 
tence, and AP is the difference between the pressure at 



ambient conditions and the pressure at coexistence. With 
AP = 1.0 X 10^ Pa, the the chemical potential change is 
A^ = 5.50 X 10~^ kBT. Note that A^ is very small. Wa- 
ter at ambient conditions is indeed close to coexistence 
with its vapor. 



B. Results 

To test whether the theory successfully addresses den- 
sity fluctuations at all length scales, a good benchmark 
is the excess chemical potential of an ideal solvophobic 
solute in a solvent as a function of its size [p] j2l]j2^ . Here 
we present results for the solvation of a hard sphere (i.e., 
spherical bubble) in water. 

The excess chemical potential A/i of a hard sphere that 
excludes solvent from a region of volume v^^ , is given by 
Eq. ( p5[ ) . We can obtain the excess chemical potential as 
a function of v^^ by sampling the size distribution of a 
"breathing" hard sphere with the Hamiltonian i?[{nfc}] 
shown in Eq. (p3|). Specifically, in a Monte Carlo tra- 
jectory for the large length scale field, each trial move 
consisted either of an attempt to flip a spin or to change 
the radius of the solute; the trial moves are accepted with 
a probability proportional to exp {—(3 AH), where AH is 
the change in [{rife}] due to the move. In order to ob- 
tain accurate statistics for all solute sizes, we have used 
umbrella sampling ||2^ . 

Figure |^ shows the excess chemical potential of a hard 
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FIG. 3. Comparison of the results of the present model 
with the predictions of the Gaussian model, Eq. (^^, and the 
results of a molecular simulation of a cavity in SPC/E water 
[14]. The horizontal line lies at the value of the surface tension 
7 of the vapor-liquid interface of water. 



sphere as a function of i?, where R is the radius of the 
spherical excluding volume. The results do exhibit some 
lattice artifacts. These artifacts, however, are surpris- 
ingly small given the fact that the cells are quite large. 
The broken lines in Fig. || reveal why the lattice artifacts 
are small. These curves show the contribution to the 
solvation free energy from the energy of the solvent, as 
given by i^L[{^^fc}], and from the energy of the "solute- 
solvent" interaction, as given by TJsbcx; {"fc}] in Eq. (23). 
Clearly, the discrete nature of our separation of length 
scales is manifest in the strong energy changes at inter- 
vals of length comparable to the grid spacing. This be- 
havior is not surprising, as the large length scale field, 
n.i, is effectively very cold (e — 6.02fcBT). The important 
point to note is that the lattice artifacts in the respective 
free-energy contributions tend to properly cancel each 
other. 

In Fig. H we compare the results of our model with the 
predictions of the Gaussian model, Eq. (p7|), and with 
the results of a molecular simulation of a hard sphere in 
SPC/E water jl^. It is seen that for small solute sizes, 
the agreement between the results of the fluid models and 
the SPC/E-simulation results is very good. The agree- 
ment is expected, since for small solutes, the large length 
scale density field remains close to its value in the un- 
perturbed fluid with [ui) « 1, as can be seen from the 
radial density profiles in Fig. ^. In this regime, our model 
reduces to the Caussian model, Eq. (|2^). Computer sim- 
ulations have shown that at small length scales, density 
fluctuations in water obey Gaussian statistics ||^. Thus, 
a Gaussian model and hence the present model success- 
fully predict the excess chemical potential of small apolar 
species in water. 

For solutes with R > 4A, the predictions of the Gaus- 
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FIG. 4. Slowly varying density {n{r)) for hard spheres 
of different size, as a function of r, where r is the distance 
to the center of the solute. This radial profile was obtained 
by averaging Ui in concentric shells of radius r and width 
Ar = O.lA. It is seen that the small solutes are in the wetting 
regime {n{r)) ~ (n(r))o ~ 1.0, whereas the larger solutes are 
in the drying regime, for which {n{r)) approaches a vapor-like 
value in the core of the solute. 



sian theory diverge from those of the full theory. The 
divergence is due to drying, as revealed in Fig. ^ The 
large length scale field (n^) approaches a vapor-like value 
in the core of the larger solutes. Gaussian models cannot 
describe this drying or depletion, because they are based 
upon a density expansion around the uniform fiuid. In 
order to describe drying, a fluid model has to support 
such a microscopic manifestation of a gas-liquid phase 
transition. 

One of the attractive features of the present model is 
that it lays bare the relative contributions to the solva- 
tion free energy. In the crossover regime, the contribu- 
tion to the Hamiltonian from the unbalancing potential, 
is very important. When we increase e' by fifty percent 
from that estimated by Eq. ( ^ ) , the results of the model 
agree very well with the simulation results over the entire 
range over which there is simulation data. In particular, 
compare Figs. ^ and ||. Evidently, orientational degrees of 
freedom result in an unbalancing potential that is larger 
than that estimated for simple fluids. With the simple 
fluid estimate of the unbalancing potential, the LCW the- 
ory ||l| overestimates the excess chemical potential in the 
crossover regime. It appears that a somewhat larger un- 
balancing potential would correct this deflciency in LCW 
theory as it does in the current model. 

It is often assumed that the excess chemical potential 
of an apolar species in water is proportional to its exposed 
surface area. Our results, as well as those of the LCW 
theory, emphasize that this is a reasonable assumption in 
the drying regime. As the solvent is near phase coexis- 
tence, the difference in chemical potential between vapor 
and liquid is small, and the work done to insert a solute 
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predominantly arises from the work to create a vapor- 
liquid interface. For small solutes, however, the excess 
chemical potential does not scale in this way. To a bet- 
ter approximation, in this regime, it is a linear function 
of fox, as can be seen from Figs. || and ||. 



VI. DISCUSSION 

We have developed a new model for a cold liquid that 
captures the effects of density fluctuations at both small 
and large length scales. This development is important, 
because many phenomena in liquids involve the interplay 
of density fluctuations in the two regimes. One exam- 
ple is capillary condensation. Another is hydrophobicity. 
Here we have focussed on the solvation of a hydropho- 
bic species. Its nature is very different at small and 
large length scales. At small length scales, solvation is 
dominated by entropic effects. In this regime, the sol- 
vent can still wet the surface of the solute, even when 
the solute is highly hydrophobic ||l|,|2l[. In contrast, at 
large length scales, solvation is dominated by energetic 
effects. In this regime, large hydrophobic objects can in- 
duce a drying transition in the solvent. Further, in the 
small length scale regime, the excess chemical potential 
scales with the volume of the solute, whereas in the large 
length scale regime, the excess chemical potential scales 
with the exposed area of the solute. The crossover be- 
havior of the solvation free energy from the wetting to 
the drying regime would seem to be of significance to 
the self-assembly of biological structures [E? 24|. In bio- 



logical systems, the size of most hydrophobic species is 
such that individual species are in the wetting regime, 
while assemblies of such species are in the drying regime. 
Water can only induce a relatively weak attraction be- 
tween two small hydrophobic species. When several of 
these species come together, however, water can induce 
a strong attraction between them. 

The crossover behavior of the solvation free energy 
also implies that the strength of the interactions between 
the hydrophobic species, depends on the configuration of 
these species. The change in the interactions manifests 
a collective effect in the solvent, and is therefore not 
simply pair decomposable. Correct simulations of self 
assembly should capture this collective effect, as could 
be done most straightforwardly with an explicit solvent 
model. While atomistic solvent models are highly lim- 
ited for this purpose, the coarse-grained model we have 
developed here should prove very useful [ pT[ . 

The effects of weak interactions have been ignored in 
this paper. Except for the movement of interfaces, that 
can be affected by small forces, weak interactions are 
not expected to induce large structural effects. Never- 
theless, their inclusion will be important for quantitative 
studies. The inclusion of weak interactions can be ac- 
complished by augmenting Eq. (^. Along with the con- 
straint factor associated with the hard core of the solute, 
the presence of a weak attractive potential, 4>(v\ between 
solute and solvent, will introduce the additional factor 
exp [— /3 J (ir0(r)] . All the analysis carried out subse- 
quent to Eq. (g) can be similarly performed in the pres- 
ence of this factor. Electrostatic interactions can also be 
incorporated, but with somewhat greater complexity. In 
this case, liquid cells (with rii — 1), must also possess a 
local polarization or dipole field, rrij. This vector field is 
Gaussian to a reasonable approximation p^ ] and there- 
fore can also be integrated out. These extensions of the 
current model are left to future work. 
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